Back to Blog | LearnableLoopAI.com | Portfolio of Projects |

versioninfo() ## Julia version

Julia Version 1.10.5
Commit 6f3fdf7b362 (2024-08-27 14:19 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 12 × Intel(R) Core(TM) i7-8700B CPU @ 3.20GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-15.0.7 (ORCJIT, skylake)
Threads: 1 default, 0 interactive, 1 GC (on 12 virtual cores)

import Pkg
Pkg.add("BenchmarkTools")

using RxInfer, BenchmarkTools, Random, LinearAlgebra, Plots, LaTeXStrings, Distributions, StableRNGs

   Resolving package versions...
  No Changes to `/workspaces/Kalman filtering and smoothing/Project.toml`
  No Changes to `/workspaces/Kalman filtering and smoothing/Manifest.toml`

Pkg.status()

Status `/workspaces/Kalman filtering and smoothing/Project.toml`
  [6e4b80f9] BenchmarkTools v1.8.0
  [31c24e10] Distributions v0.25.131
  [b964fa9f] LaTeXStrings v1.4.1
  [91a5bcdd] Plots v1.41.7
⌅ [86711068] RxInfer v3.10.1
  [860ef19b] StableRNGs v1.0.4
Info Packages marked with ⌅ have new versions available but compatibility constraints restrict them from upgrading. To see why use `status --outdated`

Kalman filtering and smoothing (Part 3)

This is an analysis of the RxInfer example at https://examples.rxinfer.com/categories/basic_examples/kalman_filtering_and_smoothing/
Some symbols have been changed
Some content has been added/modified
The preference is to make the math and code names align as much as possible
Spatial structure identifiers (e.g. vectors, matrices, cuboids)
- math
  - have a boldface e.g. \(\mathbf{x}\) or \(\boldsymbol{x}\) or \(\mathbf{A}\) or \(\mathbf{\mathbb{X}}\)
- code
  - have a underscore prefix e.g. _x or _A or _𝕏
- small case / cap case / blackboard case are generally used to discriminate between vectors, matrices, and cuboids
- start of alphabet / end of alphabet are generally used to discriminate between ‘system’ and ‘signal’
Time structure identifiers (i.e. time sequence/series)
- math
  - have a colon subscript e.g. \(x_:\) or \(x_{0:n}\) or \(x_{:n}\)
- code
  - have a ː (length mark) suffix e.g. xː
External (i.e. true, environment) states and parameters identifiers
- math
  - have a superscript * e.g. \(v^*\)
- code
  - have a superscript x e.g. vˣ
  - the ‘x’ in the code superscript is used to imitate the * superscript in the math

In the following set of examples the goal is to estimate hidden states of a Dynamical process where all hidden states are Gaussians.

We start our journey with a simple

1. multivariate Linear Gaussian State Space Model (LGSSM), which can be solved analytically. We then solve an
1. identification problem which does not have an analytical solution. Utimately, we show how RxInfer.jl can
1. deal with missing observations.

3 Handling Missing Data

An interesting case in filtering and smoothing problems is the processing of missing data. It can happen that sometimes your reading devices fail to acquire the data leading to missing observation.

Let us assume that the following model generates the data

\[\begin{aligned} {x}_t &\sim \mathcal{N}\left({x}_{t-1}, w\right) \\ {y}_t &\sim \mathcal{N}\left({x}_{t}, v \right) \end{aligned}\] where \(w = 1.0\)

with prior \({x}_0 \sim \mathcal{N}({m_{{x}_0}}, {v_{{x}_0}})\). Suppose that our measurement device fails to acquire data from time to time. In this case, instead of scalar observation \(y_t \in \mathcal{R}\) we sometimes will catch missing observations.

μᵥᵥ = 0.0
σ²ᵥᵥ = 0.1

μᵥ = 0.0
σ²ᵥ = 1.0
T = 250
missing_indices = 100:125

100:125

The Generative Process

State transition function (\(f_E\))

The state transition function provides the deterministic part of the state flow. The probabilistic part is provided by the system noise:

\[\dot{x}^*_{t} = f_{E}(t) + w = t + w\]

where

\(w \sim \mathcal{N}(\mu_W, \sigma^2_W)\)

The *s indicate that the parameters and variables are hidden and not observed.

## state transition function
function fE(;T)
    fEː = float.(collect(1:T))
    return fEː
end
myfEː = fE(T=5)

5-element Vector{Float64}:
 1.0
 2.0
 3.0
 4.0
 5.0

Observation generation function (\(g_E\))

The observation generation function provides the deterministic part of the observation. The probabilistic part is provided by the observation noise:

\[y_{t} = g_{E}(x^*_{t}) + v = \mathrm{sin}(0.05t) + v\] where

\(v \sim \mathcal{N}(\mu_V, \sigma^2_V)\)

The *s indicate that the parameters and variables are hidden and not observed.

## observation generation function
function gE(xˣː)
    gEː = sin.(0.05 * xˣː)
    return gEː
end
gE(myfEː)

5-element Vector{Float64}:
 0.04997916927067833
 0.09983341664682815
 0.14943813247359924
 0.19866933079506122
 0.24740395925452294

## Data comes from either a simulation/lab (sim|lab) OR from the field (fld)
## Data are handled either in batches (batch) OR online as individual points (point)
## Batch data accumulates either
    ## along the depth/examples dimension/axis (into the screen/page), OR
        ## typical for supervised & unsupervised learning
    ## along the time dimension/axis (down the screen page)
        ## typical for sequential decision learning (reinforcement learning & active inference)
function sim_batch_data(T; seed=123, μᵥᵥ, σ²ᵥᵥ, μᵥ, σ²ᵥ) ## simulated batch data
    rng = StableRNG(seed)
    fEː = fE(T=T)
    wː = rand(Normal(μᵥᵥ, sqrt(σ²ᵥᵥ)), T)
    xˣː = fEː + wː

    gEː = gE(xˣː)
    vː= rand(Normal(μᵥ, sqrt(σ²ᵥ)), T)
    yː = gEː + vː

    y_missingː    = similar(yː, Union{Float64, Missing}, )
    copyto!(y_missingː, yː)
    for index in missing_indices
        y_missingː[index] = missing
    end
    return fEː, xˣː, gEː, yː, y_missingː
end

sim_batch_data (generic function with 1 method)

fEː, xˣː, gEː, yː, y_missingː = sim_batch_data(
    T; 
    μᵥᵥ=μᵥᵥ, σ²ᵥᵥ=σ²ᵥᵥ, 
    μᵥ=μᵥ, σ²ᵥ=σ²ᵥ);

xˣː

250-element Vector{Float64}:
   1.3547996768376058
   1.686860549587595
   3.4270342556524307
   3.925739935761811
   4.95525999298733
   6.503616751632405
   7.038061268682471
   7.125502554706588
   9.089209701088539
   9.508955563919958
  10.90188465427425
  11.408678164919984
  12.315509492954277
   ⋮
 239.0142564013348
 239.94877013425148
 241.4968234566948
 241.95365551297226
 243.38569688085786
 244.40015255116091
 245.09073028737612
 246.31822194664505
 247.17237598747224
 247.76446243028238
 249.65131877785416
 249.78181024197582

yː

250-element Vector{Float64}:
  0.1286840362940713
 -1.0162424753368562
 -0.14128354582120928
 -0.6076284175176256
  0.44001821421305176
 -0.9056139953187743
  0.3659103660246528
 -0.28717473274632493
  0.8673813482985334
 -0.6303343105146294
  0.11119424497640196
 -1.2766482733751552
  1.4513123801338672
  ⋮
  0.4427363353645225
  0.9620157489007293
 -0.148259172600879
  0.22822378098833723
  1.4944172660785067
 -0.9855404694063241
  0.3330572010098916
  0.08496881440632997
 -0.4906913236968982
 -1.1563472928088705
  1.157323662682186
  0.11582787695490727

gEː

250-element Vector{Float64}:
  0.06768818925508203
  0.08424306389206071
  0.17051442163025998
  0.1950289793399233
  0.24523588460959617
  0.31948015150059766
  0.3446848713094134
  0.34878570137324144
  0.4389776246210882
  0.4577366072068575
  0.518498669167852
  0.539997305697964
  0.577591735842764
  ⋮
 -0.5774956663351848
 -0.538732681594585
 -0.47197481371547056
 -0.45171601443778425
 -0.3867325092065971
 -0.3394790361318947
 -0.3068047973588128
 -0.24784916487382114
 -0.206260580439247
 -0.17720668979281928
 -0.08370661341221211
 -0.077203203006163

p = plot()
p = plot(p, xˣː, gEː, label="True Signal " * L"gE_1", color=:blue)
p = scatter!(p, xˣː, yː, label="Observations", markersize=2, color=:orange)
plot(p)

The Generative Model

@model function smoothing(yː, x₀)
    v ~ Gamma(shape=0.001, scale=0.001)
    x_prior ~ Normal(mean=mean(x₀), var=var(x₀))
    local x
    xₜ₋₁ = x_prior
    for t in 1:length(yː)
        xː[t] ~ Normal(mean=xₜ₋₁, precision=1.0)
        yː[t] ~ Normal(mean=xː[t], precision=v)
        xₜ₋₁ = xː[t]
    end
end

constraints = @constraints begin
    q(x_prior, xː, yː, v) = q(x_prior, xː)q(v)q(yː)
end

Constraints: 
  q(x_prior, xː, yː, v) = q(x_prior, xː)q(v)q(yː)

x₀_prior = NormalMeanVariance(0.0, 1000.0)
initm = @initialization begin
    q(v) = Gamma(0.001, 0.001)
end

result = infer(
    model = smoothing(x₀=x₀_prior), 
    data  = (yː = y_missingː,), 
    constraints = constraints,
    initialization = initm, 
    returnvars = (xː = KeepLast(),),
    iterations = 20
);

plot(gEː, label="True signal", legend=:bottomright, color=:blue)
scatter!(
    missing_indices, 
    gEː[missing_indices], ms=2, opacity=0.75, label="Missing region", color=:blue)
plot!(
    mean.(result.posteriors[:xː]), 
    ribbon=var.(result.posteriors[:xː]), label="Estimated hidden state")